REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and 
Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person 
shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR 
FORM  TO  THE  ABOVE  ADDRESS. 


2.  REPORT  TYPE 

Journal  Article 


1.  REPORT  DATE  (DD-MM-YYYY) 
January  2013 


4.  TITLE  AND  SUBTITLE 

Thermophysical  Properties  of  Energetic  Ionic  Liquids/Nitric  Acid  Mixtures: 
Insights  from  Molecular  Dynamics  Simulations 


3.  DATES  COVERED  (From  -  To) 
January  2013-  February  2013 


5a.  CONTRACT  NUMBER 

FA9300-1  l-C-3012 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 


6.  AUTHOR(S) 


5d.  PROJECT  NUMBER 


Justin  B.  Hooper,  Grant  D.  Smith,  and  Dmitry  Bedrov 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory  (AFMC) 

AFRL/RQRP 
10  E.  Saturn  Blvd. 

Edwards  AFB  CA  93524-7680 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory  (AFMC) 

AFRL/RQR 
5  Pollux  Drive 

Edwards  AFB  CA  93524-7048 


12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Distribution  A:  Approved  for  Public  Release;  Distribution  Unlimited.  PA#13159 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 

QOPX 


8.  PERFORMING  ORGANIZATION 
REPORT  NO. 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 


11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

AFRL-RQ-ED-JA-2013-031 


13.  SUPPLEMENTARY  NOTES 

N/A 


14.  ABSTRACT 

Molecular  dynamics  (MD)  simulations  of  mixtures  of  the  room  temperature  ionic  liquids  (ILs)  1 -butyl-4  methyl  imidazolium 
[  BMIM  I'dicyanoamide  [DCA]  and  [BMIM][N03— ]  with  HN03  have  been  performed  utilizing  the  polarizable,  quantum 
chemistry  based  APPLE&P®  potential.  Experimentally  it  has  been  observed  that  [BMIM][DCA]  exhibits  hypergolic  behavior 
when  mixed  with  HN03  while  [BMIM][N03— ]  does  not.  The  structural,  thermodynamic  and  transport  properties  of  the 
IL/HN03  mixtures  have  been  determined  from  equilibrium  MD  simulations  over  the  entire  composition  range  (pure  1L  to  pure 
HN03)  based  on  bulk  simulations.  Additional  (non— equilibrium)  simulations  of  the  composition  profile  for  IL/HN03 
interfaces  as  a  function  of  time  have  been  utilized  to  estimate  the  composition  dependent  mutu  al  diffusion  coefficients  for  the 
mixtures.  The  latter  have  been  employed  in  continuum— level  simulations  in  order  to  examine  the  nature  (composition  and 
width)  of  the  IL/HN03  interfaces  on  the  millisecond  time  scale. 


16.  SECURITY  CLASSIFICATION  OF: 

a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

Unclassified 

Unclassified 

Unclassified 

17.  LIMITATION 
OF  ABSTRACT 


18.  NUMBER 
OF  PAGES 


19a.  NAME  OF 
RESPONSIBLE  PERSON 

G.  Vahjiani 


19b.  TELEPHONE  NO 

(include  area  code) 

661-525-5657 


Standard  Form 
298  (Rev.  8-98) 

Prescribed  by  ANSI 


Thermophysical  Properties  of  Energetic  Ionic  Liquids/Nitric  Acid  Mixtures: 


Insights  from  Molecular  Dynamics  Simulations 


Justin  B.  Hooper,  Grant  D.  Smith,  and  Dmitry  Bedrov 
Wasatch  Molecular,  Incorporated 
825  N.  300  W„  Salt  Lake  City,  UT  84103 


Abstract 

Molecular  dynamics  (MD]  simulations  of  mixtures  of  the  room  temperature  ionic  liquids  (ILs]  1- 
butyl-4-methyl  imidazolium  [BMIMJ/dicyanoamide  [DCA]  and  [BMIM][N03‘]  with  HNO3  have  been 
performed  utilizing  the  polarizable,  quantum  chemistry  based  APPLE&P®  potential. 
Experimentally  it  has  been  observed  that  [BMIM][DCA]  exhibits  hypergolic  behavior  when  mixed 
with  HNO3  while  [BMIM][N03'J  does  not.  The  structural,  thermodynamic  and  transport  properties 
of  the  IL/HNO3  mixtures  have  been  determined  from  equilibrium  MD  simulations  over  the  entire 
composition  range  (pure  IL  to  pure  HNO3)  based  on  bulk  simulations.  Additional  (non-equilibrium) 
simulations  of  the  composition  profile  for  IL/HNO3  interfaces  as  a  function  of  time  have  been 
utilized  to  estimate  the  composition  dependent  mutual  diffusion  coefficients  for  the  mixtures.  The 
latter  have  been  employed  in  continuum-level  simulations  in  order  to  examine  the  nature 
(composition  and  width)  of  the  IL/HNO3  interfaces  on  the  millisecond  time  scale. 
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I.  Introduction 


There  is  increasing  interest  in  use  of  room  temperature  ionic  liquids  (ILs)  as  fuels  and 
propellants.  Particularly  enticing  is  the  fact  that  RTILs  have  negligible  vapor  pressure,  making 
them  intrinsically  much  safer  than  traditional  high  vapor  pressure  and  toxic  fuels  such  as 
hydrazine.  While  most  IL  fuels  are  not  hypergolic  (i.e.,  do  not  spontaneously  combust)  when  mixed 
with  a  liquid  oxidizer  (e.g.,  nitric  acid),  there  is  a  number  of  energetic  ILs  that  do  exhibit 
hypergolicity.1  Amongst  these  are  l-butyl-4-methyl  imidazolium  [BMIM]/dicyanoamide  [DCA], 
shown  in  Figure  1.  While  [BMIM][DCA]  when  mixed  with  HNO3  is  hypergolic,  [BMIMJfNCW]  is  not.2'3 
Understanding  the  behavior  of  any  hypergolic  system  is  very  challenging  and  complicated  due  to 
the  large  manifold  of  reactions  involved  and  the  complex  interplay  between  thermophysical 
properties  and  chemical  kinetics.  When  a  fuel  droplet  comes  in  contact  with  an  oxidizer,  some 
initial  mixing  is  followed  by  the  exothermic  pre-ignition  reactions  that  begin  to  decompose  both 
fuel  and  oxidizer,  producing  various  gaseous  products.  The  pathways  for  these  reactions  involve 
highly  reactive  transient  species  that  are  difficult  to  detect  experimentally.  Changes  in  temperature 
and  species  concentrations  due  to  these  reactions  as  well  as  the  underlying  chemical  kinetics 
themselves  are  also  strongly  dependent  on  physical  factors  that  determine  heat  and  mass  transfer 
in  these  systems.  Thermophysical  properties  such  as  enthalpy  of  mixing,  species  inter-diffusion, 
and  thermal  conductivity  can  all  play  an  important  role  in  controlling  energy  dissipation  during  the 
pre-ignition  reaction  stage.2'4 

While  thermophysical  and  transport  properties  of  [BMIM][DCA]  and  [BMIMJfNCW]  have 
been  investigated  both  experimentally5'10  and  from  MD  simulations,11'18  there  are  no  data  available 
about  thermophysical  properties  of  these  liquids  mixed  with  HNO3.  Taking  into  account  that  it 
usually  takes  multiple  milliseconds  for  ignition  to  occur  in  a  hypergolic  system  there  is  a  significant 
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(from  a  molecular  scale  point  of  view]  time  period  when  thermodynamic  and  thermophysical 
processes  are  dominating  the  behavior  of  the  mixture  and  influence  the  kinetics  of  initial  reactions. 
Therefore,  we  have  carried  out  atomistic  molecular  dynamics  (MD]  simulations  of 
HN03/[BMIM][DCA]  and  HN03/[BMIM][N03']  mixtures  for  the  purposes  of  better  quantifying  the 
thermophysical  properties  of  these  mixtures  with  the  ultimate  aim  of  understanding  the  nature  of 
hypergolicity  in  ILs.  Furthermore,  modeling  of  HN03/[BMIM][DCA],  whether  at  the  continuum 
level  or  at  the  level  of  detailed  liquid-phase  reaction  pathways,  will  benefit  greatly  from  an 
improved  understanding  of  the  thermodynamics,  structure  and  transport  properties  of  the  mixture. 

II.  Force  Field  and  Simulation  Methodology 

MD  simulations  was  performed  using  the  polarizable  APPLE&P®  force  field  that  has  been 
developed  and  extensively  validated  on  a  variety  of  substances  including  ILs19  and  slightly 
modified  for  application  to  the  specific  ILs  studied  here.11  Using  this  force  field  our  previous 
simulations  showed  very  good  agreement  with  experimental  thermophysical  data  for  pure 
[BMIM][DCA]  and  [BMIM][N03‘]  RTILs.11  For  the  current  study  the  force  field  for  HNO3  was 
obtained  using  our  standard  procedure,  which  considered  all  non-bonded  interactions  and  the 
dipole  polarizability  of  all  atoms  to  be  as  transferred  from  previously  calculated  quantities  for  the 
nitrate  anion,  while  bond,  angular,  and  dihedral  interactions  were  fit  from  ab  initio  data  obtained  at 
the  MP2//aug-cc-pvDz  level.  Specific  details  about  our  fitting  procedure  can  be  found  in  Ref.  [19] 
while  the  complete  set  of  force  field  parameters  is  available  free  of  charge  upon  signing  a  user 
license  agreement. 

MD  simulations  of  bulk  [BMIM][N03‘]/HN03  and  [BMIM][DCA]/HN03  mixtures  have  been 
conducted  at  333  K  and  atmospheric  pressure  using  cubic  simulation  cell  (see  Fig.  1],  Systems  with 
a  composition  ranging  from  0.0  to  1.0  mol  fraction  IL  in  steps  of  0.2  have  been  simulated.  The  pure 
HNO3  system  contained  512  molecules  while  the  mixed  systems  employed  600,  560,  512,  or  468 
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total  molecules  in  order  of  increasing  IL  mol  fraction.  Pure  IL  data  was  taken  from  the  previous 
work,  which  employed  150  ionic  pairs.11  For  the  mixed  systems,  the  number  of  ionic  pairs 
employed  were  100,  160,  192,  and  208  for  the  mole  fractions  studied  in  increasing  IL 
concentration.  MD  simulations  were  conducted  using  the  molecular  simulation  package  Lucretius 
which  has  the  capability  to  handle  polarization  effects.  Covalent  bond  lengths  were  constrained 
using  the  velocity-Verlet  form  of  the  SHAKE  algorithm.20  The  Ewald  summation  method  was  used 
for  treatment  of  long-range  electrostatic  forces  between  partial  atomic  charges  and  between  partial 
charges  and  induced  dipoles.  A  tapering  function  was  used  to  drive  the  induced  dipole/induced 
dipole  interactions  to  zero  at  the  cutoff  of  10  A,  with  scaling  starting  at  9.3  A.  (The  form  of  the 
tapering  function  is  given  in  the  Supporting  Information)  Induced  dipoles  were  calculated  via  a 
direct  iteration  with  a  predictor  corrector  method.  A  cutoff  of  10  A  was  used  for  all  van  der  Waals 
interactions  and  the  real  part  of  electrostatic  interactions  in  the  Ewald  summation.  A  multiple  time 
step  reversible  reference  system  propagator  algorithm21  was  employed.  A  time  step  of  0.5  fs  was 
used  for  bonding,  bending,  dihedral,  and  out-of-plane  deformation  motions,  while  a  1.5  fs  time  step 
was  used  for  non-bonded  interactions  within  cutoff  radius  of  6.0  A.  Finally,  the  non-bonded 
interactions  in  the  range  between  6.0  and  10.0  A  and  reciprocal  part  of  electrostatic  interactions 
were  updated  every  3  fs.  Each  system  was  initially  equilibrated  in  the  NPT  ensemble  for  at  least  1 
ns,  while  production  runs  ranged  from  120  to  180  ns  depending  on  the  system  and  temperature. 
The  length  of  the  production  run  was  chosen  to  be  long  enough  to  allow  molecules  to  reach  a 
diffusive  regime  (when  the  ion  mean  squared  displacement  shows  a  linear  time  dependence, 
MSD[t)  ~  t),  therefore  allowing  an  accurate  estimation  of  the  self-diffusion  coefficients. 

In  addition  to  equilibrium  bulk  simulations  of  IL/HNO3  mixtures  described  above  we  also 
conducted  MD  simulations  at  non-equilibrium  conditions  in  which  we  monitor  the  process  of 
intermixing  of  IL  and  HNO3  liquids.  To  conduct  these  simulations  a  liquid-liquid  interface  has  to  be 
set  up  as  also  shown  in  Figure  1  where  RTIL  and  oxidizer  are  well  separated.  Several  adjustments 
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to  the  simulation  protocol  have  been  implemented  to  allow  these  types  of  simulations.  First,  since 
IL  and  oxidizer  are  miscible  the  preparation  of  initial  configuration  that  allows  sharp  yet  tight  (no 
vacuum]  interface  between  IL  and  oxidizer  has  to  be  prepared.  To  create  a  sharp  liquid-liquid 
interface  we  modify  parameters  for  van  der  Waals  interactions  between  IL  and  oxidizer  molecules 
to  make  them  more  repulsive.  Specifically,  we  set  dispersion  parameters  Cy  (of  the  term  -Cij/r6)  to 
zero  for  all  cross-terms  interactions  between  the  IL  and  oxidizer.  This  simple  modification  was 
enough  to  provide  a  sharp  interface  between  liquids  without  creating  any  additional  artifacts  to 
bulk  materials  away  from  the  interface.  Using  this  approach,  initial  configurations  can  be 
equilibrated  for  any  period  of  time  and  then  switched  to  simulations  with  realistic  force  field  for 
production  runs  sampling  the  intermixing  rate.  In  these  simulations,  the  simulation  cell  dimensions 
were  35  X  35  X  140  A  with  the  larger  dimension  being  perpendicular  to  the  interface.  Systems 
contained  approximately  the  same  volume  of  IL  and  oxidizer.  During  these  simulations  the 
pressure  (stress  tensor]  was  controlled  separately  in  each  direction,  allowing  the  simulation  to 
account  for  changes  in  pressure  associated  with  mixing  of  the  components.  Simulations  were 
conducted  at  333  K  and  atmospheric  pressure  control  along  the  large  dimension  (perpendicular  to 
the  interface].  Two  other  dimensions  were  fixed  at  35x35. 

III.  Simulations  of  Bulk  Mixtures 

A.  Structure 

We  began  our  analysis  of  the  HNO3/IL  mixtures  by  investigating  liquid-phase  structure  as  a 
function  of  mixture  composition.  Figures  2a-b  show  center-of-mass  BMIM-BMIM  radial 
distribution  functions  g(r]  as  a  function  of  separation  r  for  HN03/[BMIM][DCA]  and 
HN03/[BMIM][N03']  mixtures  for  the  various  compositions  investigated.  It  can  be  seen  that  for 
both  systems  there  is  a  tendency  for  non-ideal  mixing,  i.e.,  the  structure  of  the  solution  depends 
somewhat  upon  composition.  For  both  mixtures  the  position  of  the  first  peak  in  the  BMIM-BMIM 
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g(r]  moves  to  greater  separation  with  increasing  dilution  (by  HNO3].  For  HN03/[BMIM][DCA]  the 
magnitude  of  the  first  peak  also  decreases  with  increasing  dilution.  The  behavior  observed  in 
Figure  2  might  be  associated  with  strong  BMIM-HNO3  interactions,  but  Figure  3,  which  show  the 
BMIM-HNO3  radial  distribution  functions  for  the  two  mixtures  at  selected  compositions,  indicate 
that  this  is  not  the  case.  For  the  HN03/[BMIM][DCA]  mixtures  a  slight  decrease  (as  opposed  to  an 
expected  increase)  in  the  magnitude  of  the  first  peak  in  the  BMIM-HNO3  radial  distribution  function 
with  dilution  by  HNO3  can  be  seen.  Both  mixtures  show  near  ideal  mixing  behavior  with  respect  to 
BMIM-HNO3  interactions,  i.e.,  very  little  composition  dependence  is  observed  in  the  radial 
distribution  functions. 

We  have  also  investigated  the  role  of  HNO3  on  interactions  between  cations  and  anions  in 
the  mixtures  by  examining  the  cation-anion  radial  distribution  functions  as  shown  Figures  4a, b  for 
HN03/[BMIM][DCA]  and  HN03/[BMIM][N03  ],  respectively.  The  latter  shows  remarkably  little 
influence  of  the  oxidizer  on  cation-anion  correlation  in  [BMIM][N03],  while  the  former  indicates 
that  introduction  of  HNO3  promotes  stronger  nearest-neighbor  (first  radial  distribution  function 
peak)  interaction  between  BMIM  and  DCA.  Examination  of  the  anion-oxidizer  radial  distribution 
functions,  shown  in  Figures  5,  shows  very  little  dependence  of  the  nature  of  the  DCA/HNO3 
interaction  on  composition,  while  NO3'  /HNO3  interactions  become  somewhat  stronger  with 
increasing  dilution.  These  effects  are  manifested  in  DCA-DCA  pair  distribution  functions  that  are 
essentially  independent  of  composition  while  N03'/N03'  correlations  decrease  notably  with 
increasing  HNO3  content,  as  shown  in  Figures  6a-b. 

From  the  magnitude  of  the  first  peaks  in  the  radial  distribution  functions,  it  is  clear  that  the 
oxidizer  interacts  more  strongly  with  the  anions  (DCA  or  NO3')  than  with  the  BMIM  cation. 
Interestingly  the  first  peak  in  the  HNO3/NO3'  radial  distribution  function  is  high  and  narrow, 
indicating  a  strong  preference  for  specific  distance  and/or  orientations,  while  that  for  HNO3/DCA  is 
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quite  broad,  indicating  less  structure  in  this  interaction.  Atomically,  the  preferred  interaction 
between  HNO3  and  NCW  is  between  the  acidic  proton  and  the  oxygen  atoms  of  the  anion,  while  for 
HNO3  and  DCA  the  preferred  interaction  is  between  the  acidic  proton  and  the  (terminal]  nitrogen 
atoms  of  the  cyano  groups.  The  H(acidic]-O(anion]  and  H(acidic)-N(cyano]  atomic  radial 
distribution  functions  are  shown  in  Figure  7.  Interestingly,  remarkably  little  dependence  of  these 
interactions  on  the  composition  of  the  solution  is  observed.  All  radial  distribution  functions  show  a 
sharp,  high  peak  at  a  separation  of  around  two  A.  We  note  that  the  force  field  employed  in  the 
simulations  does  not  allow  for  transfer  of  the  proton  between  species,  so  the  nature  of  this 
interaction  (e.g.,  minimum  energy  position/radial  distribution  function  peak  position]  is  driven 
entirely  by  physical  (van  der  Waals  and  electrostatic]  interactions.  Integration  of  the  radial 
distribution  functions  and  multiplication  by  the  appropriate  densities  yields  the  running 
coordination  number  of  acidic  protons  near  Ofanion]  or  N(cyano]  as  shown  in  Figure  7c.  We  take  a 
separation  of  2.9  A  (see  Figures  7a-b]  as  corresponding  to  the  first  coordination  shell  of  acidic 
protons  hear  Ofanion]  or  N(cyano].  As  expected,  the  number  of  acidic  hydrogen  atoms 
coordinating  each  O(anion]  or  N(cyano]  increases  with  increasing  dilution.  For  the  Xil  =  0.2 
solutions  (80  mol%  HNO3]  the  NCh-  oxygen  atoms  and  DCA  cyano  nitrogen  atoms  are  fully 
coordinated,  i.e.,  there  is  more  than  one  on  average  of  acidic  hydrogen  atoms  in  the  first 
coordination  shell  of  these  atoms. 

B.  Density,  Volume  of  Mixing  and  Heat  of  Mixing 

The  density  and  excess  volume  for  the  two  HNO3/IL  mixtures  from  simulations  are  shown 
in  Figures  8a-b  as  a  function  of  composition.  The  density  of  the  pure  ionic  liquids  are  1.030  g/cm3 
and  1.137  g/cm3  for  [BMIM][DCA]  and  [BMIM][N03‘],  respectively.  These  values  are  in  a  good 
agreement  with  experimental  values  of  1.041  g/cm3  for  [BMIM][DCA]6  and  1.131  g/cm3  for 
[BMIM][N03']7  with  0.5%  and  1.1%  deviation,  respectively.  The  simulated  density  of  the  pure  HNO3 
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at  298  K  is  1.45  g/cm3,  compared  with  an  experimental  density22  of  1.51  g/cm3.  Mixing  of  the  ILs 
with  increasing  concentration  of  oxidizer  results  in  a  monotonic  increase  in  density.  Examination 
of  the  excess  volume  of  mixing  shows  minimal  deviation  from  ideal  mixing  behavior,  particularly 
for  the  HN03/[BMIM][N03']  system.  Figure  8c  shows  the  enthalpy  of  mixing  as  a  function  of 
composition  for  both  IL  mixtures.  Both  exhibit  exothermic  behavior,  with  the  HN03/[BMIM][N03'] 
system  showing  stronger  exothermic  character.  Combined  with  the  near-ideal  entropy  of  mixing 
(configurational  entropy)  behavior  revealed  by  the  intermolecular  radial  distribution  functions 
above,  we  can  anticipate  that  both  of  these  ILs  are  consoluble  with  HNO3.  It  is  also  clear  that  the 
temperature  of  the  systems  will  increase  upon  adiabatic  mixing,  which  is  an  important  (if  not 
crucial)  aspect  for  onset  of  ignition  in  hypergolic  mixtures.  The  extent  of  adiabatic  heating  for  these 
mixtures  is  addressed  in  Section  IV. 

C.  Self-diffusion  Coefficients 

The  rate  at  which  a  liquid  fuel  and  oxidizer  intermix  is  important  in  determining  the 
hypergolic  character  of  the  mixture.  While  intermixing  dynamics  are  investigated  in  more  detail 
below,  it  is  instructive  to  examine  the  self-diffusion  coefficients  of  the  IL  components  and  the 
oxidizer  as  a  function  of  mixture  composition  as  shown  in  Figure  9.  Figure  9  reveals  that  ion  self¬ 
diffusion  rates  are  significantly  greater  in  [BMIM][DCA]  than  in  the  [BMIM][N03‘],  consistent  with 
the  relative  viscosities  of  these  ILs  with  [BMIM][N03‘]  being  about  factor  4  more  viscous  than 
[BMIM][DCA]  at  333KA23  The  HNO3  liquid  has  a  much  lower  viscosity/greater  self-diffusion 
coefficient  than  the  ILs,  hence,  dilution  of  the  ILs  with  HNO3  results  in  a  monotonic  increase  in  ion 
self-diffusion  coefficients  and  monotonic  decrease  in  HNO3  diffusion  coefficient.  The  composition 
dependence  of  the  ion  and  HNO3  diffusion  coefficients  is  greater  in  the  HN03/[BMIM][N03']  mixture 
than  observed  in  HN03/[BMIM][DCA]  mixture.  To  first  order  this  is  consistent  with  the  greater 
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difference  in  self-diffusion  coefficients  of  HNO3  (fast)  compared  to  the  IL  (slow)  in  the  former 
system  compared  to  the  latter. 

IV.  Simulation  of  Ionic  Liquid/Oxidizer  Mixing 

In  addition  to  analysis  of  bulk  properties  under  equilibrium  conditions  discussed  above  we  have 
investigated  properties  of  mixtures  during  the  mixing  process.  Below  we  discuss  the  analysis  of 
interdiffusion  and  adiabatic  heating  processes  upon  mixing. 

A.  Diffusion 


Fig.  10a  shows  the  concentration  profile  for  HNO3  and  the  IL  the  HN03/[BMIM][DCA] 
system  after  200  ps  of  mixing  and  after  3.5  ns  of  mixing,  while  Fig.  10b  shows  the  composition 
profile  after  the  same  periods.  Similarly,  Fig.  11a  shows  the  concentration  profile  for  HNO3  and  the 
IL  the  HN03/[BMIM][N03']  system  after  200  ps  of  mixing,  after  3.0  ns  of  mixing  and  after  5.9  ns  of 
mixing,  while  Fig.  lib  shows  the  composition  profile  after  the  same  periods.  Each  profile  is  an 
average  of  three  profiles  (e.g.,  the  200  ps  profiles  are  averages  of  profiles  at  100  ps,  200  ps  and  300 
ps).  Examination  of  Figs.  10  and  11  show  obvious  asymmetry  in  mixing  at  the  interface,  indicating 
that  the  mutual  diffusion  coefficient  Dab  (A  =  oxidizer,  B  =  IL)  is  strongly  composition  dependent  for 
both  mixtures.  A  typical  diffusion  coefficient  in  these  mixtures  might  be  on  the  order  of  lxl O'9  m2/s 
(see  Fig.  9),  implying  mixing  on  a  nanometer  length  scale  on  a  nanosecond  time  scale,  consistent 
with  behavior  observed  in  Figs.  10  and  11. 


The  governing  equation  for  mutual  diffusion  (Fick's  second  law)  in  one  spatial  dimension  is 
given  as24 


dcA(r,t ) 
dt 


J_d_ 
rm  dr 


Dab{Xa) 


dcA(r,t)\ 
dr  ) 


(1) 
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(2) 


dcB{r  ,t) 
dt 


J_d_ 
rm  dr 


Dab{Xa) 


dcB(r,t ) 
dr 


X 


A 


cA(r,t ) 

cA(r,t)  +  cB(r,t) 


(3) 


where  m  =  0  for  the  slab  geometry  and  2  for  spherical  geometry.  The  dependence  of  Dab  on 
composition  is  a  priori  unknown  and  there  is  no  robust  theory  for  predicting  Dab(Xa].  However,  it 
has  been  observed  that  the  following  relationship 

l°g  Dab  =  log  Db  +  Xa  (log  Da  -  log  Db  )  (4] 

where  Da  and  Db  are  the  pure  component  self-diffusion  coefficients  works  well  in  some  materials. 
Note  that  the  self-diffusion  coefficients  of  the  species  in  the  mixtures  tend  to  follow  (at  least 
reasonably]  this  relationship  (see  Fig.  9].  Using  the  profiles  at  200  ps  as  initial  conditions,  we 
employed  eqs.  1-4  to  predict  the  concentration  and  composition  profiles  for  both  mixtures. 
Comparisons  can  be  seen  in  Figs.  10  and  11,  indicating  that  eq.  4  provides  a  reasonable  description 
of  the  composition  dependence  of  Dab-  Here  we  took  the  self-diffusion  coefficient  of  pure  HNO3 
from  Figure  9  and  used  the  average  of  the  cation  and  anion  self-diffusion  coefficients  for  the  pure  IL 
self-diffusion  coefficients,  again  taken  from  Figure  9. 

Once  the  composition  dependence  of  the  mutual  diffusion  coefficient  is  known,  it  is  possible 
to  utilize  eqs.  1-4  with  suitable  initial  and  boundary  conditions  to  prediction  the 
concentration/composition  profiles  for  any  geometry  of  interest  after  any  time  of  interest.  For 
example,  Fig.  12  shows  the  interfacial  composition  for  a  2.5  mm  radius  droplet  of  [BMIM][DCA]  and 
[BMIM][N03'],  respectively,  in  a  sea  of  HNO3  oxidizer,  after  periods  of  10  and  100  ms.  For  both 
mixtures  the  interfacial  region  is  on  the  order  of  microns  in  thickness  as  expected  given  the 
magnitude  of  the  mutual  diffusion  coefficient  (eq.  4].  It  is  notable  that  on  these  time  scales  the 
width  of  the  interfacial  region  is  not  significantly  larger  for  the  HN03/[BMIM][DCA]  mixture  than 
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for  the  HN03/[BMIM][N03']  mixture,  despite  the  fact  that  the  self-diffusion  coefficient  for 
[BMIM][DCA]  is  almost  four  times  greater  than  that  for  [BMIM][N03'].  However,  if  we  utilize  an 
imaginary  "fast"  oxidizer  with  a  self-diffusion  coefficient  four  times  that  of  HNO3,  mixing  with 
[BMIM][N03']  on  these  time  scales  occurs  much  more  rapidly,  as  shown  in  Figure  12. 

B.  Adiabatic  Heating 

Heat  generation  upon  mixing  is  another  important  property  that  can  influence  kinetics  of 
initial  reactions  and,  hence,  the  hypergolicity  of  the  mixture.  Two  approaches  were  applied  for 
calculating  heat  generation  due  to  mixing.  First,  calculations  were  based  on  bulk  mixture  heat  of 
mixing  data  (see  Section  III.B],  If  one  knows  the  heat  of  mixing  for  the  mixture  of  two  components 
(i.e.,  IL  and  oxidizer]  then  the  heat  generated  upon  their  mixing  can  be  estimated  knowing  the  final 
composition  of  the  mixture.  For  example,  for  the  [BMIM][DCA]/HN03  system  with  composition 
used  in  non-equilibrium  simulations  (Xil  =  0.23]  we  estimate  about  17  degrees  temperature 
increase.  For  [BMIM][N03]/HN03  with  similar  composition  the  corresponding  heating  is  estimated 
to  be  about  45  degrees. 

Alternatively  one  can  conduct  brute  force  simulations  of  mixing  process  in  the  NVE 
ensemble  and  directly  measure  how  much  the  system  heats  up.  While  conceptually  this  method  is 
straightforward,  special  care  must  be  taken  to  control  the  integration  uncertainty,  which  can 
otherwise  lead  to  additional  heat  generation  due  to  Hamiltonian  drift.  The  latter  can  be  measured 
by  conducting  simulations  of  a  completely  mixed  state  using  the  same  simulation  set 
up/parameters  as  in  production  run.  A  relatively  short  run  (300-500ps]  provides  a  good  estimate 
of  the  rate  of  the  Hamiltonian  drift  (if  any]  that  can  be  then  taken  into  account  (subtracted]  in  the 
analysis  of  the  production  mixing  simulation.  Using  this  approach  and  applying  correction  for  the 
Hamiltonian  drift,  a  17  degree  heating  was  obtained  from  the  non-equilibrium  simulation  of 
[BMIM][DCA]/HN03  mixing,  consistent  with  the  estimate  from  bulk  simulations  discussed  above. 
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Similarly,  the  heating  in  [BMIM]  [N03]/HN03  system  was  found  to  be  45  degrees  which  is  also 
consistent  with  the  estimate  from  bulk  simulations. 

These  simulations/estimations  of  heat  generation  due  to  physical  mixing  show  that  the 
heating  is  modest.  In  principle,  taking  into  account  that  HNO3  vaporization  temperature  is  353  K 
then  even  20-40  degrees  heating  can  be  important.  However,  it  is  important  to  realize  that  the 
relatively  narrow  intermixing  region  is  surrounded  by  the  large  amount  of  pure  IL  and  oxidizer 
phases  that  can  serve  as  heat  sinks.  In  this  case,  the  heating  of  the  mixed  layer  will  strongly  depend 
on  the  interplay  between  thermal  and  molecular  diffusivities  in  the  system.  The  latter  can  be 
characterized  by  the  Lewis  number  (Le)  and  is  the  ratio  of  the  thermal  diffusivity  of  a  fluid,  given  as 
a  =  K/(Cp*p),  to  the  molecular  diffusivity.  Taking  into  account  typical  values  for  thermal 
conductivity  (k),  heat  capacity  (Cp],  and  density  (p)  reported  in  the  literature  for  IL25  and  nitric 
acid26  we  can  estimate  a  for  our  mixture  components.  The  thermal  diffusivity  of  nitric  acid  is 
around  9xl0'8  m2/s,  while  that  for  a  number  of  various  ILs  is  in  the  range  of  6xl0'8  m2/s.  Hence, 
unlike  molecular  (mutual  or  self)  diffusivity,  we  anticipate  that  the  thermal  diffusivity  of  the  ionic 
liquid/nitric  acid  mixture  will  not  depend  strongly  on  composition,  and  will  be  on  the  order  of 
6xl0'8  m2/s.  Using  a  representative  diffusion  coefficient  of  the  IL/nitric  mixtures  of  lxlCL9  m2/s 
(which  is  on  the  high  side  of  the  range),  we  estimate  Le  >  60.  Hence,  for  viscous  fluids  such  as  ILs 
investigated  here,  the  time  scale  of  thermal  diffusion  is  much  less  than  that  of  mass  diffusion.  Even 
for  nitric  acid  itself  we  estimate  Le  ~  10. 

This  analysis  indicates  that  it  is  unlikely  that  any  noticeable  local  temperature  rise  can  be 
expected  in  the  interfacial  region  due  to  heat  of  mixing.  In  our  brute  force  mixing  simulations  we 
have  analyzed  the  local  temperature  profile  after  200-500ps.  On  these  time  scales,  IL  and  HNO3 
have  only  intermixed  in  a  relatively  small  region  of  2nm  wide.  If  the  thermal  dissipation  would  be 
much  slower  than  molecular  diffusion  then  we  should  observe  an  increase  of  temperature  in  the 
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intermixed  region  due  to  heat  of  mixing.  However,  neither  of  two  mixtures  showed  any  indication 
of  increased  temperature  in  the  intermixed  region  compared  to  the  rest  of  the  system,  indicating 
that  the  thermal  energy  generated  due  to  mixing  quickly  dissipates  to  IL  and  HNO3  phases  and 
results  in  homogeneous  heating  of  the  entire  system. 


V.  Conclusions 

Atomistic  molecular  dynamics  simulations  of  HN03/[BMIM][N03‘]  and  HN03/[BMIM][DCA] 
mixtures  using  polarizable  force  field  revealed  a  complete  consolubility  for  both  systems  in  the 
entire  range  of  concentrations.  Both  mixtures  exhibited  exothermic  behavior,  with  the 
HN03/[BMIM][N03‘]  system  showing  stronger  exothermic  character.  Analysis  of  various  pair 
correlations  functions  showed  very  little  dependence  on  mixture  composition  indicating  a  near¬ 
ideal  entropy  of  mixing.  Analysis  of  the  liquid  structure  revealed  that  the  most  prominent  difference 
between  two  mixtures  is  in  the  anion-HNCH  correlations.  The  first  peak  in  the  HNO3/NO3'  radial 
distribution  function  is  high  and  narrow,  indicating  a  strong  preference  for  specific  distance  and/or 
orientations,  while  that  for  HNO3/DCA  is  quite  broad,  indicating  less  structure  in  this  interaction. 
However,  in  both  systems  the  potential  anion  proton  acceptor  atoms  (oxygens  for  NO3'  and 
nitrogens  in  DCA]  are  well  coordinated  by  acidic  hydrogen  at  all  investigated  compositions 
indicating  that  there  no  structural  hindrances  for  initial  anion  protonation  reactions  to  occur. 

Analysis  of  molecular  diffusivity  showed  a  significant  speed  up  of  mobility  of  ions  upon 
dilution  with  nitric  acid.  We  found  that  molecular  interdiffusion  obtained  from  brute  force 
simulation  of  mixing  IL  and  oxidizer  can  be  accurately  described  by  a  simple  diffusion  equation 
utilizing  parameters  obtained  from  MD  simulations  and  simple  empirical  relation  for  calculation  of 
mutual  diffusion  coefficients.  The  developed  continuum  level  analytical  model  allowed  us  to 
estimate  concentration  profiles  on  time  and  length  scales  relevant  to  hypergolic  ignition  (multiple 
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ms  and  pm)  but  are  well  outside  of  accessibility  range  for  brute  force  atomistic  MD  simulations. 
Our  analysis  indicates  that  despite  quite  different  diffusivity/viscosity  between  two  investigated  ILs 
their  intermixing  with  nitric  acid  is  quite  similar  on  time  scales  of  10-100ms.  Finally,  the  analysis  of 
adiabatic  heating  upon  mixing  IL  and  oxidizer  revealed  that  the  released  enthalpy  of  mixing  is 
sufficient  to  allow  a  noticeable  heating  in  the  intermixed  region.  However,  due  to  intrinsically  high 
thermal  diffusivity  (compared  to  molecular  diffusivity)  of  IL  and  nitric  acid  all  heat  generated  in  the 
intermixed  layer  is  quickly  dissipating  to  non-mixed  regions.  Taking  into  account  that  on  time 
scales  10-100ms  only  a  small  fraction  of  a  fuel  droplet  has  intermixed  with  an  oxidizer  and  hence 
the  dominating  non-mixed  regions  of  IL  and  oxidizer  can  quickly  sink  the  heat  generated  in  the 
intermixed  region.  Therefore,  no  noticeable  temperature  increase  due  to  heat  of  mixing  can  be 
expected  in  the  intermixing  region. 
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Figure  Captions. 


Figure  1.  Chemical  structure  of  ions  and  nitric  acid  as  well  as  snapshots  from  equilibrium  bulk 
simulation  and  non-equilibrium  mixing  simulation  of  [BMIM][DCA]/nitric  acid  system. 

Figure  2.  BMIM-BMIM  radial  distribution  functions  for  HN03/[BMIM][DCA]  (a]  and 

HN03/[BMIM][N03‘]  (b]  mixtures. 

Figure  3.  BMIM-HNO3  radial  distribution  functions  for  HN03/[BMIM]  [DCA]  and  HN03/[BMIM][N03‘ 
]  mixtures  at  selected  compositions. 

Figure  4.  Cation-anion  radial  distribution  functions  for  HN03/[BMIM][DCA]  (a]  and 

HN03/[BMIM][N03']  [b]  mixtures. 

Figure  5.  Anion-HN03  radial  distribution  functions  for  HN03/[BMIM][DCA]  and  HN03/[BMIM][N03‘ 
]  mixtures  at  selected  compositions. 

Figure  6.  Anion-anion  radial  distribution  functions  for  HN03/[BMIM][DCA]  [a]  and 

HN03/[BMIM][N03‘]  (b]  mixtures. 

Figure  7.  The  H(acidic)-O(anion]  and  H(acidic)-N(cyano]  atomic  radial  distribution  functions  for 
HN03/[BMIM][N03']  (a]  and  HN03/[BMIM][DCA]  (b],  along  with  the  corresponding  running 
coordination  numbers  (c). 

Figure  8.  Density  (a],  excess  volume  (b]  and  heat  of  mixing  (c)  for  HN03/[BMIM][DCA]  and 
HN03/[BMIM][N03‘]  mixtures  as  a  function  of  composition. 

Figure  9.  Ion  and  oxidizer  self-diffusion  coefficients  as  a  function  of  mixture  composition  for 
HN03/[BMIM][DCA]  and  HN03/[BMIM][N03-]  mixtures. 
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Figure  10.  Time-dependent  composition  [a]  and  concentration  (b)  profiles  for  the 

HN03/[BMIM][DCA]  slab  system  after  200  ps  and  3.5  ns  of  mixing. 

Figure  11.  Time-dependent  composition  [a]  and  concentration  (b)  profiles  for  the 

HN03/[BMIM][N03']  slab  system  after  200  ps,  3.0  ns  and  5.9  ns  of  mixing. 

Figure  12.  Interfacial  profiles  for  a  2.5  mm  radius  sphere  of  [BMIM][DCA]  or  [BMIM][N03'] 
(indicated  by  arrows)  in  a  sea  of  HNO3  predicted  by  using  eqs.  1-4.  The  "fast"  oxidizer  system 
shows  mixing  of  [BMIM][N03']  with  an  oxidizer  than  has  a  self-diffusion  coefficient  four  times  that 
of  HNO3. 
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